TUTORIAL 3 : RDF

Author:Andrey V. Brukhno (andrey.brukhno{at}stfc.ac.uk)

Radial Distribution Function

Consider a system of identical particles (atoms or molecules) at constant N,V,T (canonical ensemble).

The Radial Distribution Function, or pair distribution function, g(r), describes the average density variation as a function of separation between any two particles.

To arrive at the expression for g(r), we first consider the n-particle density distribution, which is proportional to the probablity of a particular n-particle configuration (\({\mathbf r}_1, \dots, {\mathbf r}_n\)):

\[\rho^{(n)}( {\mathbf r}_1, \dots, {\mathbf r}_n ) = \frac{N!}{(N-n)!} \int \frac{\exp[-\beta U(\mathbf r^{N})]}{Z(N,V,T)} \ d{\mathbf r}_{n+1} \dots d{\mathbf r}_N\]

where the N! / (N-n)! factor accounts for all relevant permutations (due to the particles being identical); and Z(N,V,T) is the canonical partition function (normalisation constant in NVT ensemble). Setting n = 2 one can define the normalised pair distribution function

\[g({\mathbf r}_1, {\mathbf r}_2) = \frac{ \rho^{(2)}({\mathbf r}_1, {\mathbf r}_2) }{ \rho^{(1)}({\mathbf r}_1) \rho^{(1)}({\mathbf r}_2) } = \frac{V^2(N-1)}{N} \int \frac{\exp[-\beta U(\mathbf r^{N})]}{Z(N,V,T)} \ d{\mathbf r}_{3} \dots d{\mathbf r}_N\]

For an isotropic system of spherical particles one can simplify the latter formula by recognising that the radial distribution function depends only on the particle separation and should also be normalised by the area of the sphere of radius r drawn around a reference particle (“entropy on a sphere”):

\[g(r)dr = \frac{V}{4\pi r^2 N^2} \langle \sum_{i\ne j} \delta(|\Delta\mathbf r_{ij}|-r) \rangle\]

so that \(\rho(r)/\rho_b = g(r\to\infty)\to 1\), whereas \(\rho_b g(r)\) gives the local density at distance \(r\) away from a particle.

NOTE: in simulation RDF can only be determined up to \(r_{\rm max} = \frac{1}{2} \min({\rm cell\ dimensions})\)

Potential of mean force (PMF)

Potential of mean force (PMF or POMF) is defined as the work done upon reversibly bringing two particles from the state of zero-interaction into the state of full-interaction at a given distance, W(r).

PMF can be directly obtained from RDF:

\[\beta W(r) = -\ln g(r) + \ln g(\infty)\]

where the zero-interaction state corresponds to infinite separation \((r\to\infty)\).

NOTE: The PMF concept is also applicable to molecules, complexes or aggregates in solution.

Exercise 3.1 - RDF calculation

The exercises are based on TUTORIAL 2 : NPT Lennard-Jones fluid, developing it somewhat further.

Change to directory tutorial_3, where you will find the three input files for this exercise.

Below the essential amendments in CONTROL and FIELD files for this tutorial are described. Note that you don’t need to copy nor change the files from the previous tutorial!

We use reduced energy units, i.e. kT, whence Kelvin is set as energy unit in FIELD file:

UNITS K

The VDW parameters are in the units of kT and Angstrom):

VDW 1
LJ core  LJ core lj   1.0 1.0

In CONTROL file reduced temperature is set below the critical point and pressure slightly over the critical value (note that pressure is always set in katm):

temperature     1.0             # Reduced temperature T* = 1.0 < T*(CP) = 1.1876
pressure        0.02            # p* = p(katm)/0.163882576 > p*(CP)=0.1093 (reduced pressure at CP)

Collection of RDF data is invoked with the following directive in CONTROL file:

sample rdf <int> <float> <int>

where the three parameters are:

• number of grid points (bins) - set to a multiple of 10 (500)
• cutoff - set to half the (expected) minimum cell dimension (5.0)
• stride (number of MC steps) for accumulation of RDF samples - set to 1000

We also invoke collection of data for volume evolution and probability distribution, P(V), with the following line:

sample volume <int> <float>

where the two parameters are:

• stride (number of volume steps) for accumulation of samples - set to 10
• bin size - set to 10.0

Now you can re-run the NpT calculation for LJ system. Verify that RDFDAT.000 file has been created as a result. Note that RDFDAT file(s) are updated each time the stats are periodically saved, so one can check RDF data while the simulation is still being executed.

Plot the RDF in gnuplot:

[tutorial_3]$ gnuplot
gnuplot> plot 'RDFDAT.000' u 1:2 with lines

Checking the results

Questions to answer:

• Does it look reasonable?
• Any artifact observed?

You can check if the simulation converged with respect to volume (density) variations:

gnuplot> plot 'VOLDAT.000' with lines title "Volume"

You can also check the probability distribution for volume, P(V):

gnuplot> plot 'NPTDAT.000' u 1:2 with lines title "P(V)"

What might go wrong?

Questions to answer:

• Has the simulation converged w.r.t. to volume?
• Are the collected statistics sufficient (noisy data, large errors)?

You might need to perform an additional run, restarting the simulation with the latest configuration stored in REVCON file.

Before doing that, store away the current results in a new directory equil-rdf1:

[tutorial_3]$ mkdir equil-rdf1
[tutorial_3]$ cp CON* FIELD equil-rdf1/
[tutorial_3]$ mv *.0?? equil-rdf1/

Copy the obtained RECON file into CONFIG in order to restart with the latest configuration:

[tutorial_3]$ cp equil-rdf1/REVCON.000 ./CONFIG

Based on what we learnt from the equilibration run, some amendments are necessary in CONTROL file:

steps          10000000         # Number of moves to perform in simulation
equilibration      1000         # Equilibration period: statistics are gathered after this period
sample coordinate 10000         # Frequency to output configuration to trajectory file(s)
sample volume  10  2.           # Sample volume: every 10th volume move (on average); bin size
sample rdf 450 4.5 1000         # Sample RDF(s): number of bins; cutoff; stride (number of steps)

After the production run and analysis, store away the results in a new directory prod-rdf1:

[tutorial_3]$ mkdir prod-rdf1
[tutorial_3]$ cp CON* FIELD prod-rdf1/
[tutorial_3]$ mv *.0?? prod-rdf1/

Lowering density

In order to see the effect on RDF let us increase temperature (thereby, aiming at lower density). Edit CONTROL file again:

temperature     1.5            # Reduced temperature T* = 1.5 > T*(CP) = 1.1876
pressure        0.02           # p* = p(katm)/0.163882576 > p*(CP)=0.1093 (reduced pressure at CP)
...
sample rdf 400 8.0 1000
sample volume 10 10.

re-run the simulation and redo the analysis.

Again store away the current results in equil-rdf2 in case you may need to redo or restart this simulation!

Consider adding the following directives in CONTROL:

maxatmdist           5.0        # initial max atom displacement (set based on equilibration run)
acceptatmmoveupdate  1000       # stride (number of MC steps) for adjusting the max atom displacement
acceptatmmoveratio   0.4        # target acceptance ratio for atom moves [default: 0.37]

After the production run and analisys, store away the results in a new directory prod-rdf2:

[tutorial_3]$ mkdir prod-rdf2
[tutorial_3]$ cp CON* FIELD prod-rdf2/
[tutorial_3]$ mv *.0?? prod-rdf2/

Exercise 3.2 - PMF calculation

Calculate the PMF:s corresponding to the RDF:s obtained above:

gnuplot> plot [x=0.5:4.5] [y=-1:3] 'RDFDAT.000' u 1:(-log($2+0.000001)) t "PMF from RDF" w l

One can have both RDF and PMF in the same plot:

gnuplot> plot [x=0.5:4.5] [y=-1:3] 'RDFDAT.000' u 1:2 t "RDF for LJ" with lines, \
         '' u 1:(-log($2+0.000001)) t "PMF from RDF" w l

Further related tutorials

Check out TUTORIAL 6 : Free Energy (FE) and replica-exchange methods and learn how to calculate free energy differences (FED).